﻿using System;
using System.Diagnostics.Contracts;

namespace OpenLS.Spreadsheet.StandardAddIn
{
    internal class Matrix
    {
        private readonly int _columnCount;
        private readonly double[][] _data;
        private readonly int _rowCount;


        public Matrix(int rows, int columns)
        {
            _rowCount = rows;
            _columnCount = columns;
            _data = new double[rows][];
            for (int i = 0; i < rows; i++)
            {
                _data[i] = new double[columns];
            }
        }

        public Matrix(IArray array)
        {
            int columns = array.ColumnCount;
            int rows = array.RowCount;
            _rowCount = rows;
            _columnCount = columns;
            _data = new double[rows][];
            for (int i = 0; i < rows; i++)
            {
                _data[i] = new double[columns];
            }
            for (int r = 0; r < rows; r++)
            {
                for (int c = 0; c < columns; c++)
                    this[r, c] = (double)array[r, c];
            }

        }
        public Matrix(double[][] value)
        {
            _rowCount = value.Length;
            _columnCount = value[0].Length;

            for (int i = 0; i < _rowCount; i++)
            {
                if (value[i].Length != _columnCount)
                {
                    throw new ArgumentException("Argument out of range.");
                }
            }

            _data = value;
        }

        internal double[][] Array
        {
            get { return _data; }
        }

        public int RowCount
        {
            get { return _rowCount; }
        }

        public int ColumnCount
        {
            get { return _columnCount; }
        }

        public double this[int row, int column]
        {
            set { _data[row][column] = value; }

            get { return _data[row][column]; }
        }

        public Matrix Inverse
        {
            get
            {
                var m = Diagonal(_rowCount, _rowCount, 1);
                return Solve(m);
            }
        }

        public double Determinant
        {
            get { return new LuDecomposition(this).Determinant; }
        }

        public Matrix Submatrix(int[] r, int j0, int j1)
        {
            Contract.Requires(j0 <= j1);
            Contract.Requires(j0 >= 0);
            Contract.Requires(j1 >= 0);
            Contract.Assert(j0 <= _columnCount);
            Contract.Assert(j1 <= _columnCount);
            var result = new Matrix(r.Length, j1 - j0 + 1);
            var x = result.Array;
            for (var i = 0; i < r.Length; i++)
            {
                for (var j = j0; j <= j1; j++)
                {
                    if ((r[i] < 0) || (r[i] >= _rowCount))
                    {
                        throw new ArgumentException();
                    }

                    x[i][j - j0] = _data[r[i]][j];
                }
            }

            return result;
        }

        public Matrix Clone()
        {
            var result = new Matrix(_rowCount, _columnCount);
            double[][] x = result.Array;
            for (int i = 0; i < _rowCount; i++)
            {
                for (int j = 0; j < _columnCount; j++)
                {
                    x[i][j] = _data[i][j];
                }
            }

            return result;
        }


        private static Matrix Multiply(Matrix left, Matrix right)
        {
            Contract.Requires(left != null);
            Contract.Requires(right != null);
            int rows = left.RowCount;
            double[][] data = left.Array;

            if (right.RowCount != left._columnCount)
            {
                throw new ArgumentException();
            }

            int columns = right.ColumnCount;
            var result = new Matrix(rows, columns);
            double[][] x = result.Array;

            int size = left._columnCount;
            var column = new double[size];
            for (int j = 0; j < columns; j++)
            {
                for (int k = 0; k < size; k++)
                {
                    column[k] = right[k, j];
                }
                for (int i = 0; i < rows; i++)
                {
                    double[] row = data[i];
                    double s = 0;
                    for (var k = 0; k < size; k++)
                        s += row[k]*column[k];
                    x[i][j] = s;
                }
            }

            return result;
        }

        public static Matrix operator *(Matrix left, Matrix right)
        {
            Contract.Requires(left != null);
            Contract.Requires(right != null);
            return Multiply(left, right);
        }

        private Matrix Solve(Matrix m)
        {
            Contract.Requires( m != null);
            Contract.Requires( m != this);
            Contract.Requires(m.RowCount == RowCount);
            return new LuDecomposition(this).Solve(m);
        }


        private static Matrix Diagonal(int rows, int columns, double value)
        {
            Contract.Ensures(Contract.Result<Matrix>() != null);
            var result = new Matrix(rows, columns);
            var x = result.Array;
            for (int i = 0; i < Math.Min(rows, columns); i++)
                x[i][i] = value;
            return result;
        }

        public Matrix Transpose()
        {
            var result = new Matrix(_columnCount, _rowCount);
            double[][] d = result.Array;
            for (int i = 0; i < _rowCount; i++)
            {
                for (int j = 0; j < _columnCount; j++)
                {
                    d[j][i] = _data[i][j];
                }
            }

            return result;
        }
    }
}
